pacman::p_load(sf, tmap, sfdep, tidyverse, zoo, readxl, plotly)Take-home Exercise 2
Spatio-temporal Analysis of COVID-19 Vaccination Trends at the Sub-district Level, DKI Jakarta
1 Context
1.1 Problem Statement
Understanding how Covid-19 vaccination rate is distributed around DKI Jakarta and how it changes over time.
1.2 Data
| Data | Format (Type) | Description | Source |
|---|---|---|---|
| Riwayat File Vaksinasi DKI Jakarta | CSV (Aspatial) | Details the daily number of vaccinations in Jakarta by person type and sub-district (Only the first day of each month is used) | Riwatat File Vaksinasi DKI Jakarta (2022) |
| DKI Jakarta Adminstration Boundary 2019 | Shapefile (Geospatial) | DKI Jakarta administrative boundaries | Indonesia Geospasial (2019) |
Note
For the Vaksinasi DKI Jakarta data for March 2022, there was no available link to 1 March 2022 data. Thus for March 2022, the first day available, 2 March 2022, will be used instead.
2 Setting Up
2.1 Installing and Loading Packages
sf: importing and handling geospatial data
tidyverse: wrangling attribute data
sfdep: newer package of spdep (does not require conversion to sp)
tmap: prepare cartographic quality chropleth maps
zoo
readxl: read excel (xlsx) format
plotly:
2.2 DKI Jakarta Administrative Boundary 2019
2.2.1 Importing Data
To import a shapefile, use st_read() where dsn is the targeted directory and layer is the name of the files.
jkt <- st_read(
dsn = "data/geospatial",
layer = "BATAS_DESA_DESEMBER_2019_DUKCAPIL_DKI_JAKARTA"
)Reading layer `BATAS_DESA_DESEMBER_2019_DUKCAPIL_DKI_JAKARTA' from data source
`D:\bellekwang\IS415\take_home_ex\THE2\data\geospatial' using driver `ESRI Shapefile'
Simple feature collection with 269 features and 161 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 106.3831 ymin: -6.370815 xmax: 106.9728 ymax: -5.184322
Geodetic CRS: WGS 84jktis a multipolygon sf dataframe with 269 features
2.2.2 Projection
The current geodetic CRS is WGS 84 which is wrong. Using st_transform(), change the CRS from WGS 84 to DGN95 (EPSG:23845), the national projected coordinates systems of Indonesia TM-3 zone 54.1.
jkt <- st_transform(jkt, crs = 23845)
st_crs(jkt)Coordinate Reference System:
User input: EPSG:23845
wkt:
PROJCRS["DGN95 / Indonesia TM-3 zone 54.1",
BASEGEOGCRS["DGN95",
DATUM["Datum Geodesi Nasional 1995",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4755]],
CONVERSION["Indonesia TM-3 zone 54.1",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",139.5,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9999,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",200000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",1500000,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["easting (X)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["northing (Y)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Cadastre."],
AREA["Indonesia - onshore east of 138°E."],
BBOX[-9.19,138,-1.49,141.01]],
ID["EPSG",23845]]2.2.3 Data Cleaning
2.2.3.1 Outer Island
First, using tmap functions tm_shape() and tm_polygons(), plot out jkt to understand the area covered. I will also be adding the “DESA” label to identify which are the areas to exclude.
tmap_mode("view")
tm_shape(jkt) +
tm_polygons() +
tm_text("DESA")There are many outer island shown on the plot above. As seen from the plot, the areas to be excluded are: Pulau Harapan, Pulau Kelapa, Pulau Panggang, Pulau Tidung, Pulau Pari and Pulau Untung Jawa.
outer_islands <- c("PULAU HARAPAN", "PULAU KELAPA", "PULAU PANGGANG", "PULAU TIDUNG", "PULAU PARI", "PULAU UNTUNG JAWA")To remove these areas, use the filter() function.
jkt <- jkt %>%
filter(!DESA %in% outer_islands)Now just the plot mode, plot out jkt.
tmap_mode("plot")
tm_shape(jkt) +
tm_polygons()
As you can see from the plot, all the outer islands are excluded.
2.2.3.2 Unnecessary Fields
Only the first nine columns of the sf dataframe is relevant. Using the select() function, pick out the first nine columns.
jkt <- select(jkt, 1:9)2.2.4 Final jkt Data Frame
st_geometry(jkt)Geometry set for 263 features
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -3644275 ymin: 663887.8 xmax: -3606237 ymax: 701380.1
Projected CRS: DGN95 / Indonesia TM-3 zone 54.1
First 5 geometries:2.3 Riwayat File Vaksinasi DKI Jakarta
There are 12 excel files from July 2021 to June 2022 to be imported. Using read_excel() function, import the “Data Kelurahan” sheet.
july_2021 = read_excel("data/aspatial/july_2021.xlsx", sheet = 1)
aug_2021 = read_excel("data/aspatial/aug_2021.xlsx", sheet = "Data Kelurahan")
sept_2021 = read_excel("data/aspatial/sept_2021.xlsx", sheet = "Data Kelurahan")
oct_2021 = read_excel("data/aspatial/oct_2021.xlsx", sheet = "Data Kelurahan")
nov_2021 = read_excel("data/aspatial/nov_2021.xlsx", sheet = "Data Kelurahan")
dec_2021 = read_excel("data/aspatial/dec_2021.xlsx", sheet = "Data Kelurahan")
jan_2022 = read_excel("data/aspatial/jan_2022.xlsx", sheet = "Data Kelurahan")
feb_2022 = read_excel("data/aspatial/feb_2022.xlsx", sheet = "Data Kelurahan")
mar_2022 = read_excel("data/aspatial/mar_2022.xlsx", sheet = "Data Kelurahan")
apr_2022 = read_excel("data/aspatial/apr_2022.xlsx", sheet = "Data Kelurahan")
may_2022 = read_excel("data/aspatial/may_2022.xlsx", sheet = "Data Kelurahan")
june_2022 = read_excel("data/aspatial/june_2022.xlsx", sheet = "Data Kelurahan")2.3.1 Data Cleaning
I will exclude all the outer_islands identified earlier.
july_2021 <- july_2021 %>%
filter(!KELURAHAN %in% outer_islands)
aug_2021 <- aug_2021 %>%
filter(!KELURAHAN %in% outer_islands)
sept_2021 <- sept_2021 %>%
filter(!KELURAHAN %in% outer_islands)
oct_2021 <- oct_2021 %>%
filter(!KELURAHAN %in% outer_islands)
nov_2021 <- nov_2021 %>%
filter(!KELURAHAN %in% outer_islands)
dec_2021 <- dec_2021 %>%
filter(!KELURAHAN %in% outer_islands)
jan_2022 <- jan_2022 %>%
filter(!KELURAHAN %in% outer_islands)
feb_2022 <- feb_2022 %>%
filter(!KELURAHAN %in% outer_islands)
mar_2022 <- mar_2022 %>%
filter(!KELURAHAN %in% outer_islands)
apr_2022 <- apr_2022 %>%
filter(!KELURAHAN %in% outer_islands)
may_2022 <- may_2022 %>%
filter(!KELURAHAN %in% outer_islands)
june_2022 <- june_2022 %>%
filter(!KELURAHAN %in% outer_islands)All the tables now have 262 observations, 1 less than the jkt multipolygon dataframe.
2.3.1.1 Unnecessary Columns
Only the first 6 columns are important to the exercise.
july_2021 <- select(july_2021, 1:6)
aug_2021 <- select(aug_2021, 1:6)
sept_2021 <- select(sept_2021, 1:6)
oct_2021 <- select(oct_2021, 1:6)
nov_2021 <- select(nov_2021, 1:6)
dec_2021 <- select(dec_2021, 1:6)
jan_2022 <- select(jan_2022, 1:6)
feb_2022 <- select(feb_2022, 1:6)
mar_2022 <- select(mar_2022, 1:6)
apr_2022 <- select(apr_2022, 1:6)
may_2022 <- select(may_2022, 1:6)
june_2022 <- select(june_2022, 1:6)Now all the datasets only have 6 columns.
2.3.2 Combining Datasets
First, add a column to each data frame with the date using the mutate() function.
july_2021 <- july_2021 %>%
mutate(dateMonth = as.Date("2021-07-01"))
aug_2021 <- aug_2021 %>%
mutate(dateMonth = as.Date("2021-08-01"))
sept_2021 <- sept_2021 %>%
mutate(dateMonth = as.Date("2021-09-01"))
oct_2021 <- oct_2021 %>%
mutate(dateMonth = as.Date("2021-10-01"))
nov_2021 <- nov_2021 %>%
mutate(dateMonth = as.Date("2021-11-01"))
dec_2021 <- dec_2021 %>%
mutate(dateMonth = as.Date("2021-12-01"))
jan_2022 <- jan_2022 %>%
mutate(dateMonth = as.Date("2022-01-01"))
feb_2022 <- feb_2022 %>%
mutate(dateMonth = as.Date("2022-02-01"))
mar_2022 <- mar_2022 %>%
mutate(dateMonth = as.Date("2022-03-02"))
apr_2022 <- apr_2022 %>%
mutate(dateMonth = as.Date("2022-04-01"))
may_2022 <- may_2022 %>%
mutate(dateMonth = as.Date("2022-05-01"))
june_2022 <- june_2022 %>%
mutate(dateMonth = as.Date("2022-06-01"))Next, I will combine all the 12 datasets into 1 dataset using rbind() function.
vacc_data <- rbind(july_2021, aug_2021, sept_2021, oct_2021, nov_2021, dec_2021, jan_2022, feb_2022, mar_2022, apr_2022, may_2022, june_2022)2.3.3 Final Dataset
head(vacc_data)# A tibble: 6 × 7
`KODE KELURAHAN` `WILAYAH KOTA` KECAMATAN KELUR…¹ SASARAN BELUM…² dateMonth
<chr> <chr> <chr> <chr> <dbl> <dbl> <date>
1 <NA> <NA> <NA> TOTAL 7739060 5041111 2021-07-01
2 3172051003 JAKARTA UTARA PADEMANGAN ANCOL 20393 13272 2021-07-01
3 3173041007 JAKARTA BARAT TAMBORA ANGKE 25785 16477 2021-07-01
4 3175041005 JAKARTA TIMUR KRAMAT JATI BALE K… 25158 18849 2021-07-01
5 3175031003 JAKARTA TIMUR JATINEGARA BALI M… 8683 5743 2021-07-01
6 3175101006 JAKARTA TIMUR CIPAYUNG BAMBU … 22768 15407 2021-07-01
# … with abbreviated variable names ¹KELURAHAN, ²`BELUM VAKSIN`3 Choropleth Mapping and Analysis
3.1 Monthly Vaccination Rates at Sub-district Level
3.1.1 Computing Monthly Vaccination Rates
SASARAN (Target - target people to be vaccinated)
BELUM VAKSIN (Yet to be vaccinated)
To compute the monthly vaccination rate,
\[ Monthly Vaccination Rate = \frac{(SASARAN - BELUMVAKSIN)}{SUM(SASARAN, BELUMVAKSIN)} \]
vacc_data <- vacc_data %>%
mutate(monthly_vacc_rate = (SASARAN - `BELUM VAKSIN`)/(SASARAN + `BELUM VAKSIN`))3.1.2 Combining Dataset and Shapefile
vacc_july_2021 <- vacc_data %>%
filter(dateMonth == "2021-07-01")
jkt_july_2021 <- left_join(jkt, vacc_july_2021, by = c("DESA" = "KELURAHAN" ))3.2 Choropleth Mapping
Using the qtm() function, plot out the choropleth map for July 2021.
qtm(jkt_july_2021, "monthly_vacc_rate")
As seen above, there are missing monthly vaccinated rate values for some sub-districts. There is a mismatch and lack of vaccination data for these areas.
I will leave them as missing values instead of giving a monthly vaccination rate of 0 to avoid confusion with actual areas with recorded 0 rates.
3.2.1 Other Months
Repeat the above steps for the next 11 months.
Code
#AUGUST 2021
vacc_aug_2021 <- vacc_data %>%
filter(dateMonth == "2021-08-01")
jkt_aug_2021 <- left_join(jkt, vacc_aug_2021, by = c("DESA" = "KELURAHAN" ))
#SEPTEMBER 2021
vacc_sept_2021 <- vacc_data %>%
filter(dateMonth == "2021-09-01")
jkt_sept_2021 <- left_join(jkt, vacc_sept_2021, by = c("DESA" = "KELURAHAN" ))
#OCTOBER 2021
vacc_oct_2021 <- vacc_data %>%
filter(dateMonth == "2021-10-01")
jkt_oct_2021 <- left_join(jkt, vacc_oct_2021, by = c("DESA" = "KELURAHAN" ))
#NOVEMBER 2021
vacc_nov_2021 <- vacc_data %>%
filter(dateMonth == "2021-11-01")
jkt_nov_2021 <- left_join(jkt, vacc_nov_2021, by = c("DESA" = "KELURAHAN" ))
#DECEMBER 2021
vacc_dec_2021 <- vacc_data %>%
filter(dateMonth == "2021-12-01")
jkt_dec_2021 <- left_join(jkt, vacc_dec_2021, by = c("DESA" = "KELURAHAN" ))
#JANUARY 2022
vacc_jan_2022 <- vacc_data %>%
filter(dateMonth == "2022-01-01")
jkt_jan_2022 <- left_join(jkt, vacc_jan_2022, by = c("DESA" = "KELURAHAN" ))
#FEBRUARY 2022
vacc_feb_2022 <- vacc_data %>%
filter(dateMonth == "2022-02-01")
jkt_feb_2022 <- left_join(jkt, vacc_feb_2022, by = c("DESA" = "KELURAHAN" ))
#MARCH 2022
vacc_mar_2022 <- vacc_data %>%
filter(dateMonth == "2022-03-02")
jkt_mar_2022 <- left_join(jkt, vacc_mar_2022, by = c("DESA" = "KELURAHAN" ))
#APRIL 2022
vacc_apr_2022 <- vacc_data %>%
filter(dateMonth == "2022-04-01")
jkt_apr_2022 <- left_join(jkt, vacc_apr_2022, by = c("DESA" = "KELURAHAN" ))
#MAY 2022
vacc_may_2022 <- vacc_data %>%
filter(dateMonth == "2022-05-01")
jkt_may_2022 <- left_join(jkt, vacc_may_2022, by = c("DESA" = "KELURAHAN" ))
#JUNE 2022
vacc_june_2022 <- vacc_data %>%
filter(dateMonth == "2022-06-01")
jkt_june_2022 <- left_join(jkt, vacc_june_2022, by = c("DESA" = "KELURAHAN" ))Code
qtm(jkt_july_2021, "monthly_vacc_rate", title = "July 2021")
Code
qtm(jkt_aug_2021, "monthly_vacc_rate", title = "Aug 2021")
Code
qtm(jkt_sept_2021, "monthly_vacc_rate", title = "Sept 2021")
Code
qtm(jkt_oct_2021, "monthly_vacc_rate", title = "Oct 2021")
Code
qtm(jkt_nov_2021, "monthly_vacc_rate", title = "Nov 2021")
Code
qtm(jkt_dec_2021, "monthly_vacc_rate", title = "Dec 2021")
Code
qtm(jkt_jan_2022, "monthly_vacc_rate", title = "Jan 2022")
Code
qtm(jkt_feb_2022, "monthly_vacc_rate", title = "Feb 2022")
Code
qtm(jkt_mar_2022, "monthly_vacc_rate", title = "Mar 2022")
Code
qtm(jkt_apr_2022, "monthly_vacc_rate", title = "Apr 2022")
Code
qtm(jkt_may_2022, "monthly_vacc_rate", title = "May 2022")
Code
qtm(jkt_june_2022, "monthly_vacc_rate", title = "June 2022")
3.3 Spatial Patterns
4 Local Gi* Analysis
4.1 Computing Local Gi* of Monthly Vaccination Rate
To compute local Gi* statistics, it requires the neighbour list (using st_contiguity())and the wts (using st_inverse_distance()).
Using jkt_july_2021 as an example:
jkt_july_2021 = jkt_july_2021 %>% drop_na(monthly_vacc_rate)
jkt_july_2021_idw <- jkt_july_2021 %>%
mutate(
nb = st_contiguity(geometry),
wts = st_inverse_distance(
nb, geometry, scale = 1, alpha = 1
),
.before = 1
)HCSA_july_2021 <- jkt_july_2021_idw %>%
mutate(local_Gi = local_gstar_perm(
monthly_vacc_rate, nb, wt, nsim = 99),
.before = 1) %>%
unnest(local_Gi)
HCSA_july_2021Simple feature collection with 252 features and 26 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -3644275 ymin: 663887.8 xmax: -3606237 ymax: 701380.1
Projected CRS: DGN95 / Indonesia TM-3 zone 54.1
# A tibble: 252 × 27
gi_star e_gi var_gi p_value p_sim p_fol…¹ skewn…² kurtosis nb wts
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <nb> <lis>
1 1.76 0.00401 1.11e-7 0.0778 0.1 0.05 0.319 -0.231 <int> <dbl>
2 3.02 0.00394 8.95e-8 0.00251 0.04 0.02 0.709 1.30 <int> <dbl>
3 -0.570 0.00399 1.08e-7 0.569 0.62 0.31 0.481 0.148 <int> <dbl>
4 -1.13 0.00399 8.53e-8 0.260 0.24 0.12 0.0785 0.0614 <int> <dbl>
5 2.40 0.00394 6.52e-8 0.0163 0.06 0.03 0.166 0.501 <int> <dbl>
6 1.05 0.00398 1.14e-7 0.293 0.28 0.14 0.177 0.660 <int> <dbl>
7 0.933 0.00394 8.73e-8 0.351 0.38 0.19 0.00187 -0.00350 <int> <dbl>
8 -0.0584 0.00398 1.04e-7 0.953 0.9 0.45 0.303 -0.740 <int> <dbl>
9 0.607 0.00402 7.31e-8 0.544 0.56 0.28 0.141 -0.365 <int> <dbl>
10 -0.138 0.00392 7.94e-8 0.890 0.96 0.48 0.0895 -0.476 <int> <dbl>
# … with 242 more rows, 17 more variables: OBJECT_ID <dbl>, KODE_DESA <chr>,
# DESA <chr>, KODE <dbl>, PROVINSI <chr>, KAB_KOTA <chr>, KECAMATAN.x <chr>,
# DESA_KELUR <chr>, JUMLAH_PEN <dbl>, `KODE KELURAHAN` <chr>,
# `WILAYAH KOTA` <chr>, KECAMATAN.y <chr>, SASARAN <dbl>,
# `BELUM VAKSIN` <dbl>, dateMonth <date>, monthly_vacc_rate <dbl>,
# geometry <MULTIPOLYGON [m]>, and abbreviated variable names ¹p_folded_sim,
# ²skewness4.2 Displaying Gi* Maps
(only significant values: p<0.05)
tmap_mode("plot")
map1 <- tm_shape(HCSA_july_2021) +
tm_fill("gi_star", palette = "-RdBu") +
tm_borders(alpha = 0.5) +
tm_view(set.zoom.limits = c(6,8)) +
tm_layout(main.title = "Gi* of July 2021",
main.title.size = 0.8)
map2 <- tm_shape(HCSA_july_2021) +
tm_fill("p_value",
breaks = c(0, 0.001, 0.01, 0.05, 1),
labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") +
tm_borders(alpha = 0.5) +
tm_layout(main.title = "p-value of Gi*",
main.title.size = 0.8)
tmap_arrange(map1, map2, ncol = 2)
4.3 Other Months
Repeat the previous steps to the other 11 months.
4.3.1 Aug 2021
Code
jkt_aug_2021 = jkt_aug_2021 %>% drop_na(monthly_vacc_rate)
jkt_aug_2021_idw <- jkt_aug_2021 %>%
mutate(
nb = st_contiguity(geometry),
wts = st_inverse_distance(
nb, geometry, scale = 1, alpha = 1
),
.before = 1
)
HCSA_aug_2021 <- jkt_aug_2021_idw %>%
mutate(local_Gi = local_gstar_perm(
monthly_vacc_rate, nb, wt, nsim = 99),
.before = 1) %>%
unnest(local_Gi)
tmap_mode("plot")
map1 <- tm_shape(HCSA_aug_2021) +
tm_fill("gi_star", palette = "-RdBu") +
tm_borders(alpha = 0.5) +
tm_view(set.zoom.limits = c(6,8)) +
tm_layout(main.title = "Gi* of August 2021",
main.title.size = 0.8)
map2 <- tm_shape(HCSA_aug_2021) +
tm_fill("p_value",
breaks = c(0, 0.001, 0.01, 0.05, 1),
labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") +
tm_borders(alpha = 0.5) +
tm_layout(main.title = "p-value of Gi*",
main.title.size = 0.8)
tmap_arrange(map1, map2, ncol = 2)
4.3.2 September 2021
Code
jkt_sept_2021 = jkt_sept_2021 %>% drop_na(monthly_vacc_rate)
jkt_sept_2021_idw <- jkt_sept_2021 %>%
mutate(
nb = st_contiguity(geometry),
wts = st_inverse_distance(
nb, geometry, scale = 1, alpha = 1
),
.before = 1
)
HCSA_sept_2021 <- jkt_sept_2021_idw %>%
mutate(local_Gi = local_gstar_perm(
monthly_vacc_rate, nb, wt, nsim = 99),
.before = 1) %>%
unnest(local_Gi)
tmap_mode("plot")
map1 <- tm_shape(HCSA_sept_2021) +
tm_fill("gi_star", palette = "-RdBu") +
tm_borders(alpha = 0.5) +
tm_view(set.zoom.limits = c(6,8)) +
tm_layout(main.title = "Gi* of September 2021",
main.title.size = 0.8)
map2 <- tm_shape(HCSA_sept_2021) +
tm_fill("p_value",
breaks = c(0, 0.001, 0.01, 0.05, 1),
labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") +
tm_borders(alpha = 0.5) +
tm_layout(main.title = "p-value of Gi*",
main.title.size = 0.8)
tmap_arrange(map1, map2, ncol = 2)
4.3.3 October 2021
Code
jkt_oct_2021 = jkt_oct_2021 %>% drop_na(monthly_vacc_rate)
jkt_oct_2021_idw <- jkt_oct_2021 %>%
mutate(
nb = st_contiguity(geometry),
wts = st_inverse_distance(
nb, geometry, scale = 1, alpha = 1
),
.before = 1
)
HCSA_oct_2021 <- jkt_oct_2021_idw %>%
mutate(local_Gi = local_gstar_perm(
monthly_vacc_rate, nb, wt, nsim = 99),
.before = 1) %>%
unnest(local_Gi)
tmap_mode("plot")
map1 <- tm_shape(HCSA_oct_2021) +
tm_fill("gi_star", palette = "-RdBu") +
tm_borders(alpha = 0.5) +
tm_view(set.zoom.limits = c(6,8)) +
tm_layout(main.title = "Gi* of October 2021",
main.title.size = 0.8)
map2 <- tm_shape(HCSA_oct_2021) +
tm_fill("p_value",
breaks = c(0, 0.001, 0.01, 0.05, 1),
labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") +
tm_borders(alpha = 0.5) +
tm_layout(main.title = "p-value of Gi*",
main.title.size = 0.8)
tmap_arrange(map1, map2, ncol = 2)
4.3.4 November 2022
Code
jkt_nov_2021 = jkt_nov_2021 %>% drop_na(monthly_vacc_rate)
jkt_nov_2021_idw <- jkt_nov_2021 %>%
mutate(
nb = st_contiguity(geometry),
wts = st_inverse_distance(
nb, geometry, scale = 1, alpha = 1
),
.before = 1
)
HCSA_nov_2021 <- jkt_nov_2021_idw %>%
mutate(local_Gi = local_gstar_perm(
monthly_vacc_rate, nb, wt, nsim = 99),
.before = 1) %>%
unnest(local_Gi)
tmap_mode("plot")
map1 <- tm_shape(HCSA_nov_2021) +
tm_fill("gi_star", palette = "-RdBu") +
tm_borders(alpha = 0.5) +
tm_view(set.zoom.limits = c(6,8)) +
tm_layout(main.title = "Gi* of November 2021",
main.title.size = 0.8)
map2 <- tm_shape(HCSA_nov_2021) +
tm_fill("p_value",
breaks = c(0, 0.001, 0.01, 0.05, 1),
labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") +
tm_borders(alpha = 0.5) +
tm_layout(main.title = "p-value of Gi*",
main.title.size = 0.8)
tmap_arrange(map1, map2, ncol = 2)
4.3.5 December 2021
Code
jkt_dec_2021 = jkt_dec_2021 %>% drop_na(monthly_vacc_rate)
jkt_dec_2021_idw <- jkt_dec_2021 %>%
mutate(
nb = st_contiguity(geometry),
wts = st_inverse_distance(
nb, geometry, scale = 1, alpha = 1
),
.before = 1
)
HCSA_dec_2021 <- jkt_dec_2021_idw %>%
mutate(local_Gi = local_gstar_perm(
monthly_vacc_rate, nb, wt, nsim = 99),
.before = 1) %>%
unnest(local_Gi)
tmap_mode("plot")
map1 <- tm_shape(HCSA_dec_2021) +
tm_fill("gi_star", palette = "-RdBu") +
tm_borders(alpha = 0.5) +
tm_view(set.zoom.limits = c(6,8)) +
tm_layout(main.title = "Gi* of December 2021",
main.title.size = 0.8)
map2 <- tm_shape(HCSA_dec_2021) +
tm_fill("p_value",
breaks = c(0, 0.001, 0.01, 0.05, 1),
labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") +
tm_borders(alpha = 0.5) +
tm_layout(main.title = "p-value of Gi*",
main.title.size = 0.8)
tmap_arrange(map1, map2, ncol = 2)
4.3.6 January 2022
Code
jkt_jan_2022 = jkt_jan_2022 %>% drop_na(monthly_vacc_rate)
jkt_jan_2022_idw <- jkt_jan_2022 %>%
mutate(
nb = st_contiguity(geometry),
wts = st_inverse_distance(
nb, geometry, scale = 1, alpha = 1
),
.before = 1
)
HCSA_jan_2022 <- jkt_jan_2022_idw %>%
mutate(local_Gi = local_gstar_perm(
monthly_vacc_rate, nb, wt, nsim = 99),
.before = 1) %>%
unnest(local_Gi)
tmap_mode("plot")
map1 <- tm_shape(HCSA_jan_2022) +
tm_fill("gi_star", palette = "-RdBu") +
tm_borders(alpha = 0.5) +
tm_view(set.zoom.limits = c(6,8)) +
tm_layout(main.title = "Gi* of January 2022",
main.title.size = 0.8)
map2 <- tm_shape(HCSA_jan_2022) +
tm_fill("p_value",
breaks = c(0, 0.001, 0.01, 0.05, 1),
labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") +
tm_borders(alpha = 0.5) +
tm_layout(main.title = "p-value of Gi*",
main.title.size = 0.8)
tmap_arrange(map1, map2, ncol = 2)
4.3.7 February 2022
Code
jkt_feb_2022 = jkt_feb_2022 %>% drop_na(monthly_vacc_rate)
jkt_feb_2022_idw <- jkt_feb_2022 %>%
mutate(
nb = st_contiguity(geometry),
wts = st_inverse_distance(
nb, geometry, scale = 1, alpha = 1
),
.before = 1
)
HCSA_feb_2022 <- jkt_feb_2022_idw %>%
mutate(local_Gi = local_gstar_perm(
monthly_vacc_rate, nb, wt, nsim = 99),
.before = 1) %>%
unnest(local_Gi)
tmap_mode("plot")
map1 <- tm_shape(HCSA_feb_2022) +
tm_fill("gi_star", palette = "-RdBu") +
tm_borders(alpha = 0.5) +
tm_view(set.zoom.limits = c(6,8)) +
tm_layout(main.title = "Gi* of February 2022",
main.title.size = 0.8)
map2 <- tm_shape(HCSA_feb_2022) +
tm_fill("p_value",
breaks = c(0, 0.001, 0.01, 0.05, 1),
labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") +
tm_borders(alpha = 0.5) +
tm_layout(main.title = "p-value of Gi*",
main.title.size = 0.8)
tmap_arrange(map1, map2, ncol = 2)
4.3.8 March 2022
Code
jkt_mar_2022 = jkt_mar_2022 %>% drop_na(monthly_vacc_rate)
jkt_mar_2022_idw <- jkt_mar_2022 %>%
mutate(
nb = st_contiguity(geometry),
wts = st_inverse_distance(
nb, geometry, scale = 1, alpha = 1
),
.before = 1
)
HCSA_mar_2022 <- jkt_mar_2022_idw %>%
mutate(local_Gi = local_gstar_perm(
monthly_vacc_rate, nb, wt, nsim = 99),
.before = 1) %>%
unnest(local_Gi)
tmap_mode("plot")
map1 <- tm_shape(HCSA_mar_2022) +
tm_fill("gi_star", palette = "-RdBu") +
tm_borders(alpha = 0.5) +
tm_view(set.zoom.limits = c(6,8)) +
tm_layout(main.title = "Gi* of March 2022",
main.title.size = 0.8)
map2 <- tm_shape(HCSA_mar_2022) +
tm_fill("p_value",
breaks = c(0, 0.001, 0.01, 0.05, 1),
labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") +
tm_borders(alpha = 0.5) +
tm_layout(main.title = "p-value of Gi*",
main.title.size = 0.8)
tmap_arrange(map1, map2, ncol = 2)
4.3.9 April 2022
Code
jkt_apr_2022 = jkt_apr_2022 %>% drop_na(monthly_vacc_rate)
jkt_apr_2022_idw <- jkt_apr_2022 %>%
mutate(
nb = st_contiguity(geometry),
wts = st_inverse_distance(
nb, geometry, scale = 1, alpha = 1
),
.before = 1
)
HCSA_apr_2022 <- jkt_apr_2022_idw %>%
mutate(local_Gi = local_gstar_perm(
monthly_vacc_rate, nb, wt, nsim = 99),
.before = 1) %>%
unnest(local_Gi)
tmap_mode("plot")
map1 <- tm_shape(HCSA_apr_2022) +
tm_fill("gi_star", palette = "-RdBu") +
tm_borders(alpha = 0.5) +
tm_view(set.zoom.limits = c(6,8)) +
tm_layout(main.title = "Gi* of April 2022",
main.title.size = 0.8)
map2 <- tm_shape(HCSA_apr_2022) +
tm_fill("p_value",
breaks = c(0, 0.001, 0.01, 0.05, 1),
labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") +
tm_borders(alpha = 0.5) +
tm_layout(main.title = "p-value of Gi*",
main.title.size = 0.8)
tmap_arrange(map1, map2, ncol = 2)
4.3.10 May 2022
Code
jkt_may_2022 = jkt_may_2022 %>% drop_na(monthly_vacc_rate)
jkt_may_2022_idw <- jkt_may_2022 %>%
mutate(
nb = st_contiguity(geometry),
wts = st_inverse_distance(
nb, geometry, scale = 1, alpha = 1
),
.before = 1
)
HCSA_may_2022 <- jkt_may_2022_idw %>%
mutate(local_Gi = local_gstar_perm(
monthly_vacc_rate, nb, wt, nsim = 99),
.before = 1) %>%
unnest(local_Gi)
tmap_mode("plot")
map1 <- tm_shape(HCSA_may_2022) +
tm_fill("gi_star", palette = "-RdBu") +
tm_borders(alpha = 0.5) +
tm_view(set.zoom.limits = c(6,8)) +
tm_layout(main.title = "Gi* of May 2022",
main.title.size = 0.8)
map2 <- tm_shape(HCSA_may_2022) +
tm_fill("p_value",
breaks = c(0, 0.001, 0.01, 0.05, 1),
labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") +
tm_borders(alpha = 0.5) +
tm_layout(main.title = "p-value of Gi*",
main.title.size = 0.8)
tmap_arrange(map1, map2, ncol = 2)
4.3.11 June 2022
Code
jkt_june_2022 = jkt_june_2022 %>% drop_na(monthly_vacc_rate)
jkt_june_2022_idw <- jkt_june_2022 %>%
mutate(
nb = st_contiguity(geometry),
wts = st_inverse_distance(
nb, geometry, scale = 1, alpha = 1
),
.before = 1
)
HCSA_june_2022 <- jkt_june_2022_idw %>%
mutate(local_Gi = local_gstar_perm(
monthly_vacc_rate, nb, wt, nsim = 99),
.before = 1) %>%
unnest(local_Gi)
tmap_mode("plot")
map1 <- tm_shape(HCSA_june_2022) +
tm_fill("gi_star", palette = "-RdBu") +
tm_borders(alpha = 0.5) +
tm_view(set.zoom.limits = c(6,8)) +
tm_layout(main.title = "Gi* of June 2022",
main.title.size = 0.8)
map2 <- tm_shape(HCSA_june_2022) +
tm_fill("p_value",
breaks = c(0, 0.001, 0.01, 0.05, 1),
labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") +
tm_borders(alpha = 0.5) +
tm_layout(main.title = "p-value of Gi*",
main.title.size = 0.8)
tmap_arrange(map1, map2, ncol = 2)
4.4 Statistical Conclusion
5 Emerging Hot Spot Analysis
5.1 Mann-Kendall Test
(using spatio-temporal local Gi* values)
5.1.1 Creating Time Series Cube
vacc_data = vacc_data %>% drop_na(c("KODE KELURAHAN", "WILAYAH KOTA", "KECAMATAN"))
vacc_data_st <- spacetime(vacc_data, jkt, ".DESA", "dateMonth")is_spacetime_cube(vacc_data_st)[1] FALSE5.1.2 Computing Spatial Weights
5.2 Temporal Trends
(select 3 sub-districts)
5.3 EHSA Map
(only significant values)